
# =============================================================
# File: App_D_sector.R
# Purpose: Produces results from Appendix Section D: Additional Sector Fixed Effects Models
# Paper: Foreignness as an Asset: European Carbon Regulation and the Relocation Threat among Multinational Firms
# Author: Patrick Bayer, patrick.bayer@strath.ac.uk
# Date: 5 October 2022
#
# Data: ./data.csv
#
# Technical disclaimer:
# All analyses in R version 4.1.2 (2021-11-01)
# RStudio 2022.07.1 Build 554 ("Spotted Wakerobin" Release (7872775e, 2022-07-22) for Windows)
# Windows 10 Enterprise, 64-bit
# 12th Gen Intel(R) Core(TM) i7-1270P 2.20 GHz with 32GB RAM
# =============================================================

library(lmtest)
library(sandwich)
library(MASS)
library(MatchIt)
library(cem)
library(cobalt)
library(cowplot)
library(ggpubr)
library(marginaleffects)
library(AER)
library(tidyverse)

# Load data
df <- read_csv(file = "./data.csv")

# Trim sample down to true within firm sample
# (1) Limit sample to European MNCs: no allocation data for domestically-owned operations for non-EU MNCs
# (2) Exclude MNCs with only foreign-owned operations as there is no within MNC variation

df$sample <- ifelse(df$mnc.eu==1 & df$foreign.share<1,1,0)
df <- df[df$sample==1,]
df$foreign <- as.factor(df$foreign) # set variable as categorical

# =============================================================
# APPENDIX D:  Additional sector fixed effects models
# =============================================================

# In the main text, sectors are defined in terms of activity codes as used in the EU ETS
# This is useful because some countries used permit allocation rules along these sectoral lines
# Sectors can alternatively be defined according to NACE codes, the EU's official industry standard classification system
# Here, I replicate models with sector fixed effects for NACE codes from 1-4 digit levels

# Note: Model (1) does not use sector fixed effects, hence I replicate results for Models (2) and (3) from main text only

# Code various NACE classification levels
# Source: https://ec.europa.eu/competition/mergers/cases/index/nace_all.html

# Create different sets of NACE code variables
df$NACE2 <- substr(df$NACE,1,nchar(df$NACE)-3) 
df$NACE2 <- ifelse(df$NACE2=="",NA,df$NACE2)

# Create NACE3 variable
df$NACE3 <- substr(df$NACE,1,nchar(df$NACE)-1)   
df$NACE3 <- ifelse(df$NACE3=="",NA,df$NACE3)

# Create NACE4 variable
df$NACE4 <- ifelse(df$NACE=="",NA,df$NACE)

# Create NACE1 variable
df$NACE2num <- as.numeric(df$NACE2)

df$NACE1 <- ifelse(df$NACE2num<=3, "A",
ifelse(df$NACE2num>3 & df$NACE2num<=9, "B",
ifelse(df$NACE2num>9 & df$NACE2num<=33, "C",
ifelse(df$NACE2num>33 & df$NACE2num<=35, "D",
ifelse(df$NACE2num>35 & df$NACE2num<=39, "E",
ifelse(df$NACE2num>39 & df$NACE2num<=43, "F",
ifelse(df$NACE2num>43 & df$NACE2num<=47, "G",
ifelse(df$NACE2num>47 & df$NACE2num<=53, "H",
ifelse(df$NACE2num>53 & df$NACE2num<=56, "I",
ifelse(df$NACE2num>56 & df$NACE2num<=63, "J",
ifelse(df$NACE2num>63 & df$NACE2num<=66, "K",
ifelse(df$NACE2num>66 & df$NACE2num<=68, "L",
ifelse(df$NACE2num>68 & df$NACE2num<=75, "M",
ifelse(df$NACE2num>75 & df$NACE2num<=82, "N",
ifelse(df$NACE2num>82 & df$NACE2num<=84, "O",
ifelse(df$NACE2num>84 & df$NACE2num<=85, "P",
ifelse(df$NACE2num>85 & df$NACE2num<=88, "Q",
ifelse(df$NACE2num>88 & df$NACE2num<=93, "R",
ifelse(df$NACE2num>93 & df$NACE2num<=96, "S",
ifelse(df$NACE2num>96 & df$NACE2num<=98, "T",
ifelse(df$NACE2num>98 & df$NACE2num<=99, "U", ""
)))))))))))))))))))))

# =============================================================
# Model (2) with different NACE-level sector FEs
# =============================================================

# Model(2) for NACE1
m2 <- lm(logDV~foreign+logAlloc0+logEmit0+as.factor(iso2)+as.factor(NACE1)+as.factor(BVD), data=df)
summary(m2)

# Observations
nobs(m2)
length(m2$xlevels$`as.factor(NACE1)`)
length(m2$xlevels$`as.factor(iso2)`)
length(m2$xlevels$`as.factor(BVD)`)

# Marginal effects
(exp(summary(m2)$coef[["foreign1","Estimate"]])-1)*100

# Standard errors
(exp(coefci(m2, vcov. = vcovHC, type="HC0")["foreign1",])-1)*100 # robust SEs


# Model(2) for NACE2
m2 <- lm(logDV~foreign+logAlloc0+logEmit0+as.factor(iso2)+as.factor(NACE2)+as.factor(BVD), data=df)
summary(m2)

# Observations
nobs(m2)
length(m2$xlevels$`as.factor(NACE2)`)
length(m2$xlevels$`as.factor(iso2)`)
length(m2$xlevels$`as.factor(BVD)`)

# Marginal effects
(exp(summary(m2)$coef[["foreign1","Estimate"]])-1)*100

# Standard errors
(exp(coefci(m2, vcov. = vcovHC, type="HC0")["foreign1",])-1)*100 # robust SEs


# Model (2) for NACE3
m2 <- lm(logDV~foreign+logAlloc0+logEmit0+as.factor(iso2)+as.factor(NACE3)+as.factor(BVD), data=df)
summary(m2)

# Observations
nobs(m2)
length(m2$xlevels$`as.factor(NACE3)`)
length(m2$xlevels$`as.factor(iso2)`)
length(m2$xlevels$`as.factor(BVD)`)

# Marginal effects
(exp(summary(m2)$coef[["foreign1","Estimate"]])-1)*100

# Standard errors
(exp(coefci(m2, vcov. = vcovHC, type="HC0")["foreign1",])-1)*100 # robust SEs


# Model (2) for NACE4
m2 <- lm(logDV~foreign+logAlloc0+logEmit0+as.factor(iso2)+as.factor(NACE4)+as.factor(BVD), data=df)
summary(m2)

# Observations
nobs(m2)
length(m2$xlevels$`as.factor(NACE4)`)
length(m2$xlevels$`as.factor(iso2)`)
length(m2$xlevels$`as.factor(BVD)`)

# Marginal effects
(exp(summary(m2)$coef[["foreign1","Estimate"]])-1)*100

# Standard errors
(exp(coefci(m2, vcov. = vcovHC, type="HC0")["foreign1",])-1)*100 # robust SEs


# =============================================================
# Model (3) with different NACE-level sector FEs
# =============================================================

# Exact matching on 4-digit NACE code and firm
df$NACE <- ifelse(df$NACE=="",NA,df$NACE) # Recode empty strings as NA 

m <- matchit(foreign ~ NACE + BVD, data=df[is.na(df$NACE)==FALSE,], method="exact")
df.match <- match.data(m)

# Model (3) for NACE1
m3 <- lm(logDV~foreign+logAlloc0+logEmit0+as.factor(iso2)+as.factor(NACE1)+as.factor(BVD), data=df.match, weights=weights)
summary(m3)

# Observations
nobs(m3)
length(m3$xlevels$`as.factor(NACE1)`)
length(m3$xlevels$`as.factor(iso2)`)
length(m3$xlevels$`as.factor(BVD)`)

# Marginal effects
(exp(summary(m3)$coef[["foreign1","Estimate"]])-1)*100

# Standard errors
(exp(coefci(m3, vcov. = vcovHC, type="HC0")["foreign1",])-1)*100 # robust SEs


# Model (3) for NACE2
m3 <- lm(logDV~foreign+logAlloc0+logEmit0+as.factor(iso2)+as.factor(NACE2)+as.factor(BVD), data=df.match, weights=weights)
summary(m3)

# Observations
nobs(m3)
length(m3$xlevels$`as.factor(NACE2)`)
length(m3$xlevels$`as.factor(iso2)`)
length(m3$xlevels$`as.factor(BVD)`)

# Marginal effects
(exp(summary(m3)$coef[["foreign1","Estimate"]])-1)*100

# Standard errors
(exp(coefci(m3, vcov. = vcovHC, type="HC0")["foreign1",])-1)*100 # robust SEs


# Model (3) for NACE3
m3 <- lm(logDV~foreign+logAlloc0+logEmit0+as.factor(iso2)+as.factor(NACE3)+as.factor(BVD), data=df.match, weights=weights)
summary(m3)

# Observations
nobs(m3)
length(m3$xlevels$`as.factor(NACE3)`)
length(m3$xlevels$`as.factor(iso2)`)
length(m3$xlevels$`as.factor(BVD)`)

# Marginal effects
(exp(summary(m3)$coef[["foreign1","Estimate"]])-1)*100

# Standard errors
(exp(coefci(m3, vcov. = vcovHC, type="HC0")["foreign1",])-1)*100 # robust SEs


# Model (3) for NACE4
m3 <- lm(logDV~foreign+logAlloc0+logEmit0+as.factor(iso2)+as.factor(NACE4)+as.factor(BVD), data=df.match, weights=weights)
summary(m3)

# Observations
nobs(m3)
length(m3$xlevels$`as.factor(NACE4)`)
length(m3$xlevels$`as.factor(iso2)`)
length(m3$xlevels$`as.factor(BVD)`)

# Marginal effects
(exp(summary(m3)$coef[["foreign1","Estimate"]])-1)*100

# Standard errors
(exp(coefci(m3, vcov. = vcovHC, type="HC0")["foreign1",])-1)*100 # robust SEs


# =============================================================
# Models with country-sector FEs
# =============================================================

# Create country-sector FE variables
df$ctrysector1 <- ifelse(is.na(df$NACE1)==FALSE,paste(df$iso2,df$NACE1), NA)
df$ctrysector2 <- ifelse(is.na(df$NACE2)==FALSE,paste(df$iso2,df$NACE2), NA)
df$ctrysector3 <- ifelse(is.na(df$NACE3)==FALSE,paste(df$iso2,df$NACE3), NA)
df$ctrysector4 <- ifelse(is.na(df$NACE4)==FALSE,paste(df$iso2,df$NACE4), NA)


# =============================================================
# Model (2) with country-sector FEs at different NACE code hierarchies
# =============================================================

# Model(2) with NACE1 for country-sector FEs
m2 <- lm(logDV~foreign+logAlloc0+logEmit0+as.factor(ctrysector1)+as.factor(BVD), data=df)
summary(m2)

# Observations
nobs(m2)
length(m2$xlevels$`as.factor(ctrysector1)`)
length(m2$xlevels$`as.factor(BVD)`)

# Marginal effects
(exp(summary(m2)$coef[["foreign1","Estimate"]])-1)*100

# Standard errors
(exp(coefci(m2, vcov. = vcovHC, type="HC0")["foreign1",])-1)*100 # robust SEs



# Model(2) with NACE2 for country-sector FEs
m2 <- lm(logDV~foreign+logAlloc0+logEmit0+as.factor(ctrysector2)+as.factor(BVD), data=df)
summary(m2)

# Observations
nobs(m2)
length(m2$xlevels$`as.factor(ctrysector2)`)
length(m2$xlevels$`as.factor(BVD)`)

# Marginal effects
(exp(summary(m2)$coef[["foreign1","Estimate"]])-1)*100

# Standard errors
(exp(coefci(m2, vcov. = vcovHC, type="HC0")["foreign1",])-1)*100 # robust SEs



# Model(2) with NACE3 for country-sector FEs
m2 <- lm(logDV~foreign+logAlloc0+logEmit0+as.factor(ctrysector3)+as.factor(BVD), data=df)
summary(m2)

# Observations
nobs(m2)
length(m2$xlevels$`as.factor(ctrysector3)`)
length(m2$xlevels$`as.factor(BVD)`)

# Marginal effects
(exp(summary(m2)$coef[["foreign1","Estimate"]])-1)*100

# Standard errors
(exp(coefci(m2, vcov. = vcovHC, type="HC0")["foreign1",])-1)*100 # robust SEs




# Model(2) with NACE4 for country-sector FEs
m2 <- lm(logDV~foreign+logAlloc0+logEmit0+as.factor(ctrysector4)+as.factor(BVD), data=df)
summary(m2)

# Observations
nobs(m2)
length(m2$xlevels$`as.factor(ctrysector4)`)
length(m2$xlevels$`as.factor(BVD)`)

# Marginal effects
(exp(summary(m2)$coef[["foreign1","Estimate"]])-1)*100

# Standard errors
(exp(coefci(m2, vcov. = vcovHC, type="HC0")["foreign1",])-1)*100 # robust SEs




# =============================================================
# Model (3) with country-sector FEs at different NACE code hierarchies
# =============================================================


# Exact matching on 4-digit NACE code and firm
df$NACE <- ifelse(df$NACE=="",NA,df$NACE) # Recode empty strings as NA 

m <- matchit(foreign ~ NACE + BVD, data=df[is.na(df$NACE)==FALSE,], method="exact")
df.match <- match.data(m)


# Model(3) with NACE1 for country-sector FEs
m3 <- lm(logDV~foreign+logAlloc0+logEmit0+as.factor(ctrysector1)+as.factor(BVD), data=df.match, weights=weights)
summary(m3)

# Observations
nobs(m3)
length(m3$xlevels$`as.factor(ctrysector1)`)
length(m3$xlevels$`as.factor(BVD)`)

# Marginal effects
(exp(summary(m3)$coef[["foreign1","Estimate"]])-1)*100

# Standard errors
(exp(coefci(m3, vcov. = vcovHC, type="HC0")["foreign1",])-1)*100 # robust SEs



# Model(3) with NACE2 for country-sector FEs
m3 <- lm(logDV~foreign+logAlloc0+logEmit0+as.factor(ctrysector2)+as.factor(BVD), data=df.match, weights=weights)
summary(m3)

# Observations
nobs(m3)
length(m3$xlevels$`as.factor(ctrysector2)`)
length(m3$xlevels$`as.factor(BVD)`)

# Marginal effects
(exp(summary(m3)$coef[["foreign1","Estimate"]])-1)*100

# Standard errors
(exp(coefci(m3, vcov. = vcovHC, type="HC0")["foreign1",])-1)*100 # robust SEs
(exp(coefci(m3, vcov. = vcovHC, type="HC0", level=.9)["foreign1",])-1)*100 # robust SEs




# Model(3) with NACE3 for country-sector FEs
m3 <- lm(logDV~foreign+logAlloc0+logEmit0+as.factor(ctrysector3)+as.factor(BVD), data=df.match, weights=weights)
summary(m3)

# Observations
nobs(m3)
length(m3$xlevels$`as.factor(ctrysector3)`)
length(m3$xlevels$`as.factor(BVD)`)

# Marginal effects
(exp(summary(m3)$coef[["foreign1","Estimate"]])-1)*100

# Standard errors
(exp(coefci(m3, vcov. = vcovHC, type="HC0")["foreign1",])-1)*100 # robust SEs
(exp(coefci(m3, vcov. = vcovHC, type="HC0", level=.9)["foreign1",])-1)*100 # robust SEs



# Model(3) with NACE4 for country-sector FEs
m3 <- lm(logDV~foreign+logAlloc0+logEmit0+as.factor(ctrysector4)+as.factor(BVD), data=df.match, weights=weights)
summary(m3)

# Observations
nobs(m3)
length(m3$xlevels$`as.factor(ctrysector4)`)
length(m3$xlevels$`as.factor(BVD)`)

# Marginal effects
(exp(summary(m3)$coef[["foreign1","Estimate"]])-1)*100

# Standard errors
(exp(coefci(m3, vcov. = vcovHC, type="HC0")["foreign1",])-1)*100 # robust SEs


# =============================================================
#                       END OF CODE
# =============================================================